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Abstract 

Meshfree finite difference methods for the Poisson equation approx- 
imate the Laplace operator on a point cloud. Desirable are positive 
stencils, i.e. all neighbor entries are of the same sign. Classical least 
squares approaches yield large stencils that are in general not positive. 
We present an approach that yields stencils of minimal size, which 
are positive. We provide conditions on the point cloud geometry, so 
that positive stencils always exist. The new discretization method is 
compared to least squares approaches in terms of accuracy and com- 
putational performance. 

1 Introduction 

The numerical approximation of the Poisson equation is a fundamental task 
encountered in many applications. Often it appears as a subproblem in a 
more complex computation, for instance as a projection step in the simu- 
lation of incompressible flows [TJ. Finite difference methods approximate 
the equation on a finite number of points. If the points can be placed on 
a regular grid, the approximation is simple and yields symmetric matrices. 
However, in many cases a regular point distribution is not possible or de- 
sired. Examples are the explicit representation of complex geometries, or the 
point positions may be given by the application, for instance from scattered 
measurements or in particle methods [8] . If the points are distributed irreg- 
ularly, neighborhood relations have to be established. This could be done by 
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constructing a mesh. However, meshing can be costly, and thus may not be 
desired in applications with time-dependent geometries. Instead, meshfree 
neighborhood criteria can be defined, and meshfree finite difference stencils 
constructed. Consistency conditions for stencils are derived in Sect. [2j We 
are interested in approaches that yield M-matrices, as explained in Sect. 
This requires the stencils to be positive. Least squares approaches, outlined 
in Sect.Sl in general fail to yield positive stencils. In Sect.[5]we present a new 
approach, based on sign constrained linear minimization, that yields posi- 
tive stencils. In Sect. [6] we provide conditions on the point cloud geometry, 
so that positive stencils are guaranteed to exist. Further conditions, derived 
in Sect. ensure an M-matrix structure. In Sect. [8] the new approach is 
compared to classical methods by numerical experiments. 

2 Meshfree Finite Differences for the Poisson Equa- 
tion 

Consider the Poisson equation to be solved inside a domain Q C K d 



where Tjj U Tjy = Oil. Let a point cloud X = {xi, . . . , x n } C fl be given, 
which consists of interior points Ij C and boundary points Xj, C dfl. 
The point cloud is meshfree, i.e. no information about connection of points 
is provided. Meshfree finite difference approaches convert problem ([I]) into 
a linear system 



where the vector u contains approximations to the values u(xi). The i-th 
row of the matrix A consists of the stencil corresponding to the point a?j. 
We assume that ([TJ admits a unique smooth solution. 

2.1 Consistent Derivative Approximation 

Consider a function u G C 2 (Q C R rf ,R). We wish to approximate Au(xq) 
using the function values of a finite number of points in a circular neighbor- 
hood (xq,xi, . . . ,x m ) G B(xo,r), where B(xo,r) = {x G Q : \\x — xq\\ < r}. 
Define the distance vectors Xi = x^ — xq Vz = 0, . . . , m. The function value 
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at each neighboring point u(xi) can be expressed by a Taylor expansion 

u(xi) = u(xq) + Vu(x Q ) ■ Xi + \V 2 u(xq) : (a3j ■ xf) + . 

We use the matrix scalar product A : B = ^ - AyBij. The error in the 
expansion is of order a = 0(r 3 ). A linear combination with coefficients 
(s ,... ,s m ) equals 



^ Si«(xi) = u(a3o) E s » + V-u(a;o) • ^ s, 



i=0 



vi=0 



vi=l 
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This approximates the Laplacian, i.e. s i u ( x i) = Au(a; ) + 0(r 3 ), if 

exactness for constant, linear and quadratic functions is satisfied 



i=0 







^ a^Sj = , ^ (a?j • xj) Si = 21 . 
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Definition 1. A stencil (sq, . . . , s m ) to a set of points (xq, x\, . . . , x m ) £ 
B(xQ,r) is called consistent (with the Laplace operator), if the constraints 
(|3|) are satisfied. 

The linear and quadratic constraints in ([3]) can be formulated as a linear 
system of equations 

V ■ s = b , (4) 

where V £ R^™ is the Vandermonde matrix given by Xi,...,x m , and 
s € W 71 is the stencil vector. In 2d, with X{ = (xi,y~i), the system reads as 
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The constant constraint 



yields sq = — YaLi s «- Neumann boundary points can be treated in a similar 
manner. Approximating -§^(xq) by Yli=o s i u ( x i) leads to the constraints 
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For each point, a meshfree finite difference approximation consists of two 
steps: First, define which points are its neighbors. Typically, more neigh- 
bors than constraints are chosen. Second, select a stencil. If (J3j) is underde- 
termined, a minimization problem is formulated to select a unique stencil. 

A set of neighbors around a central point is called in general configura- 
tion, if the Vandermonde matrix V has full rank. For m neighboring points, 
one has no solution if m < k, infinitely many solutions if m > k, and one 
solution if m = k. In this case, the stencil s = V^ 1 ■ b can be computed by 
elimination, or by formulas for the determinant of a multivariate Vander- 
monde matrix [12] . If the points are not in general configuration, e.g. for 
regular grids, the above rules may fail. A solution can exist for m < k 
(e.g. 5-point stencil in 2d), and no solution may exist for m > k (see exam- 
ple in [151 p. 59]). The concept of general configuration is unhandy and too 
strict. Solutions may be acceptable, even if V does not have full rank. The 
geometric condition presented in Sect. [6] ensures the existence of stencils. 

2.2 Minimal and Positive Stencils 

Definition 2. A consistent stencil (so, • • • , s m ) is called minimal, if m < k. 

Minimal stencils are beneficial for the sparsity of the system matrix, 
resulting in a lower memory consumption and a faster solution of system 
(|2J) • The total number of neighboring points is proportional to the effort of 
applying the matrix to a vector, which is proportional to the time for one 
step of a (semi-)iterative linear solver. Of course, the few neighbors have to 
be chosen wisely, to preserve good convergence rates of iterative solvers. The 
results in [TBI E] indicate that this is the case with the presented approach. 

Remark 1. For minimal stencils, it is impossible that the stencil values de- 
pend continuously on the point positions. Consider six points around a 
central point (in 2d), five of which are selected neighbors. Consider a con- 
tinuous movement of one of the neighbors and the sixth point, such that at 
the end these two points have swapped their positions, without them ever 
being in the same place. At some instance during this movement, the sixth 
point has to become a neighbor, resulting in a jump in the stencil values. 

Remark 2. If used in a particle method, the lack of smoothness in minimal 
stencils (Rem. [JJ may lead to a non-conservative scheme. In an isolated 
Poisson solver and in particle methods that are not conservative by con- 
struction (such as the finite pointset method [8]) the advantage of optimal 
sparsity often outweighs this drawback. 
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Definition 3. A consistent stencil (sq, . . . , s m ) is called positive, if 
si, ■ ■ ■ ,s m > 0. Due to ([3]) and ([6]), this implies for the central point so < 0. 

Positive stencils yields the system matrix in ([2]) to be an L-matrix 
(Def. H]), which gives rise to an M-matrix structure (see Sect. ED- The de- 
sirability of positive stencils has been pointed out by Demkowicz, Karafiat 
and Liszka [2j. Classical approaches do in general not yield positive stencils. 
An "optimal star selection" [3] makes positive stencils likely, but they are 
not guaranteed (see Fig. [5|). Fiirst and Sonar derive topological conditions 
on point clouds for positive least squares stencils in Id [6]. In Sect. [5] we 
present a strategy that approximates the Poisson equation Jl| on a point 
cloud by minimal positive stencils. In Sect. [6] conditions on a point cloud are 
presented (in 2d and 3d) which guarantee the existence of positive stencils. 

3 M-Matrices 

Meshfree finite difference matrices are in general non-symmetric. Consider 
two points Xi and Xj, each being a neighbor of the other, and a third point 
Xk which is a neighbor of Xi but not a neighbor of Xj. Since each stencil 
entry depends on all its neighbors, the point x^ influences the matrix entry 
a«, but not the matrix entry aji. 

The negative Laplace operator in (pQ) is positive definite. For non- 
symmetric matrices, we have to ask for slightly less than positive definite- 
ness. A property which implies a maximum principle and the convergence 
of linear solvers, is the M-matrix structure. 

Definition 4. A square matrix A = (aij)ij € M nxn is called Z-matrix, if 
<Hj < Vi ^ j. A Z-matrix is called L-matrix, if an > Vi. 

We write A > for a^ > Vi,j. The same notation applies to vectors. 

Definition 5. A regular matrix A is called inverse positive, if A^ 1 > 0. 

Definition 6. A Z-matrix is called M-matrix, if it is inverse positive. 

We use the M-matrix property, since it yields a sufficient condition for in- 
verse positivity. There are inverse positive matrices that are not M-matrices, 
so another approach would be to employ alternative characterizations of in- 
verse positive matrices pjj. 
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3.1 Benefits of an M-Matrix Structure 



The Poisson equation satisfies maximum principles. For instance, consider 
(PQ) with Dirichlet boundary conditions only. If / < and g < 0, then the 
solution satisfies u < [I]. A discretization by an M- matrix mimics this 
property in a discrete maximum principle. 

Theorem 1. Let A be an M-matrix. Then Ax < implies x < 0. Con- 
versely, a Z-matrix satisfying Ax < =4> x < is an M-matrix. 

Proof. A is an M-matrix, thus A -1 > by definition. Let y = Ax. Then 
x = A~ 1 y. The component- wise inequalities A^ 1 > and y < imply 
x < 0. The reverse statement is proved in |13|. p. 29]. □ 

Theorem 2. // A is an M-matrix, and D its diagonal part, then p(I — 
D~ l A) < 1, thus the Jacobi and the Gaufi-Seidel iteration converge. 

Proof. The convergence of the Jacobi iteration is given in [7]. The GauB- 
Seidel convergence follows from the Stein-Rosenberg-Theorem [20]. □ 

The performance of multigrid methods for meshfree finite difference ma- 
trices has been investigated in [16\ [T7] . Further benefits of an M-matrix 
structure with respect to linear solvers can be found in [20j . 

3.2 A Sufficient Condition for an M-Matrix Structure 

Conditions that imply the M-matrix property are required, since the inverse 
matrix is typically not directly available. Let the unknowns be labeled by 
an index set I. We consider square matrices A G R /x/ . 

Definition 7. The graph G(A) of a matrix A is defined by G{A) = {(i,j) £ 
I x / : aij ^ 0}. The index i £ I is called connected to j G I, if a chain i = 
io,i\, . . . , ik~i,ik = j 6 I exists, such that (i u -i,i u ) G G{A) Vi/ = 1, . . . , k. 

Remark 3. For a finite difference matrix, each index i £ I corresponds to a 
point X{. The index i £ I being connected to j £ I means that the point Xi 
connects (indirectly) to the point Xj via stencil entries. 

Definition 8. A finite difference matrix is called essentially irreducible if 
every point is connected to a Dirichlet boundary point. 

Remark 4. A finite difference matrix that is not essentially irreducible, is 
singular, since the points that are not connected to a Dirichlet point form a 
singular submatrix. 
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Definition 9. A matrix A G ~R IxI is called essentially diagonally dominant, 
if it is weakly diagonally dominant (Vi El: \au\ > and every 

point i G / is connected to a point j G J which satisfies the strict diagonal 
dominance relation \djj\ > J2k^j \ a jk\- 

Theorem 3. An L-matrix arising as a finite difference discretization of ([1]) 
is essentially diagonally dominant, if it is essentially irreducible. 

Proof. For an L-matrix the constant relation in ([3]) implies the weak diag- 
onal dominance relation for every interior and Neumann point. Each row 
corresponding to a Dirichlet point satisfies the strict diagonal dominance 
relation. □ 

Theorem 4. An essentially diagonally dominant L-matrix is an M-matrix. 

Proof. The proof is given in (TJ p. 153]. □ 

If problem ([1]) can be discretized by positive stencils and every point is 
connected to a Dirichlet point, then the resulting matrix is an M-matrix. 



4 Least Squares Approaches 

Classical approaches for meshfree derivative approximation are moving least 
squares methods, based on scattered data interpolation [9J, and local approx- 
imation methods, based on generalized finite difference methods [11]. Their 
application to meshfree settings has been analyzed in j3j 110j . Differences 
between moving and local approaches have been investigated in [15J. 

Around a central point xq, points inside a radius r are considered. A 
distance weight function w(5) is defined, which is small for 5 > r. We 
consider interpolating approaches with w(5) = 5~ a . Each neighboring 
point Xi is assigned a weight Wi = w (\\xi — £Co H2)- A unique stencil is 
defined via a quadratic minimization problem 

n 2 

minV^-, s.t. V ■ s = b . (7) 

i=l 

Using W = diag(wj), its solution is 

s = WV T (VWV T )~ 1 b . (8) 

1 If lima^o w(S) exists, the approach is called approximating, otherwise interpolating. 
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After the weights Wi are evaluated, the k x k matrix VWV T has to be 
set up, the linear system (VWV T ) ■ v = b to be solved, and the product 
s = WV T ■ v to be computed. This requires k(k + l)m + ^ floating point 
operations [15j P- 150]. Least squares approaches do not yield minimal 
stencils, unless exactly k neighbors are considered. In general, they also do 
not yield positive stencils. 

Example 1. Consider xq = (0,0) and 6 neighbors on the unit circle Xi = 
(cos(f^),sin(f^)), where fa, . . . ,<p 6 ) = (0,1,2,3,0.1,0.2) (see Fig. p. 
Since all neighbors have the same distance from xq, the distance weight 
function does not play a role. Formula (|8|) yields the non-positive least 
squares stencil s = (0.846,1.005,0.998,1.003,0.312,-0.164). However, the 
configuration admits a positive stencil, namely s = (1, 1, 1, 1, 0, 0). 



5 Linear Minimization Approach 

Least squares approaches do not guarantee positive stencils. As motivated 
in Sect. 12.21 we wish to allow positive stencils only. Hence, we enforce 
positivity, i.e. we search for solutions in the polyhedron 

P = {s G R m : V ■ s = b, s > 0} . (9) 

This is the feasibility problem of linear optimization [19]. In Sect. [6] we show 
under which conditions solutions exist. If P is nonvoid and not degenerate, 
there are infinitely many feasible stencils. To single out a unique stencil we 
formulate a linear minimization problem 

m 

minV— , s.t. V-s = b, s>0, (10) 

*H Wi 

i=i 

where the weights Wi = w(\\xi — x^\\) are defined by an appropriately de- 
caying (see Thm. [6]) non-negative distance weight function w. Problem (|10p 
is a linear program (LP) in standard form. It is bounded, since we have 
imposed sign constraints and the weights Wi are all non-negative. 

Theorem 5. If the polyhedron is nonvoid, then the linear minimization 
approach (jlQI) yields minimal positive stencils. 

Proof. The sign constraints in (|10p ensure that the selected stencil is positive. 
The existence of a minimal solution is ensured by the fundamental theorem 
of linear programming [19]. If the LP (|10p has a solution, then it also has 
a basic solution, in which at most k of the m stencil entries Si are different 
from zero. □ 
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Figure 1: Non-positive Figure 2: Minimal positive stencil for various val- 
LSQ stencil ues of a 

In contrast to least squares methods, the LP approach yields nonzero 
values only for a few selected points, and a continuous dependence on the 
point positions is not possible (see Rem. [[]). 

Remark 5. One could ask why not remain with a least squares problem, and 
additionally impose sign constraints. While this would be a valid approach 
(the solution is obtained by Karush-Kuhn- Tucker methods [E]), it would 
have the worst of both worlds. The solution would not depend continu- 
ously on the point cloud geometry whenever the sign constraints are active, 
and the resulting stencil would not be minimal. When sign constraints are 
imposed, linear minimization is preferable. 

5.1 Solving the Linear Programs 

For every interior point consider a set of candidate points (m > k). A 
basic solution of (I10p is computed. Only the nonzero stencil values enter 
the Poisson matrix. We refer to this approach as minimal positive stencil 
(MPS) method. The LPs (jlOh are small, but they have to be solved for every 
interior point. To our knowledge, there are no general results about efficient 
methods for such small LPs, especially considering the special structure of 
the Vandermonde matrix. A numerical comparison of various methods has 
been presented in [15, p. 148]. Simplex methods perform best for the arising 
LPs. A basis change corresponds to one stencil point replacing another. 
The theoretical worst case performance of simplex methods is not observed. 
Typical runs find the solution in about 1.5k steps, resulting in a complexity 
of 0(k 2 m), which equals the effort of least squares approaches (see Sect. dj. 
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5.2 Geometric Interpretation of Minimal Positive Stencils 

The MPS method forms a compromise between selecting neighbors close by 
and distributed nicely (see Def. 1X0]) around the central point. How much 
preference is given to which objective depends on the locality parameter a 
in the distance weight function w(5) = \5\~ a . 

Theorem 6. The MPS method (|10p only leads to reasonable results, if the 
distance weight function decays faster than \d~\~ 2 ■ 

Proof. Summing over the diagonal in the quadratic constraints in ([3]) yields 
the relation Y^iLi ll^i|li s i = 2d. If w(5) decays faster than \5\~ 2 , points 
close to the central point are given preference. If w(S) = \b~\~ 2 , the LP 
(fTOj) is degenerate. If w(5) decays slower than \5\~ 2 , the approach selects 
points far away from the central point, possibly resulting in "checkerboard" 
instabilities. □ 

The dependence of the MPS stencil on a is shown in Fig. [2^1 From the 
candidate points in the circle, five neighbors are selected. While for a = 1 
far away points are selected, a G {3, 5} yields nearby points. For a = 3 
smaller angles are more important, for a = 5 smaller distances. Note that 
the MPS method never selects neighbors which are not distributed around 
the central point (as defined in Sect. [6]), even if those are the k closest points. 

Remark 6. For regular grids, the MPS method selects standard finite dif- 
ference stencils. For instance, for a regular Cartesian grid, the standard 
5-point (2d), respectively 7-point (3d) stencils are obtained. In these cases, 
the basic solution is degenerate, i.e. some of the basis variables are zero. 

5.3 Neighborhood Criteria 

The circular neighborhood criterion yields a large number of neighbors (un- 
less r is very small). Also, it does not guarantee positive stencils, as the 
example in Fig. [3] shows. The selected neighbors are marked bold. Neigh- 
bors with non-positive stencil values are indicated by a white center. The 
presented methods can also be based on other neighborhood criteria. 

Defining a neighborhood via the a Delaunay triangulation (tetrahedriza- 
tion in 3d) [18] . yields significantly fewer neighbors. However, the construc- 
tion is a meshing procedure, hence it is often undesirable in a meshfree 

2 The figures are 2d for simplicity of presentation. The MPS method applies directly 
to, and shows its strength, in 3d. 
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Figure 3: Circu- Figure 4: Delau- Figure 5: Four Figure 6: MPS 
lar neighborhood nay neighbors quadrants neighborhood 

contextH Fig. 0] shows that a Delaunay neighborhood constructed from a 
Voronoi tessellation may yield non-positive stencils. In addition, it is pos- 
sible that not enough neighbors are selected. If, in the example, the two 
negative points were removed, only four neighbors would be defined. 

The four quadrant criterion [3] (eight sectors in 3d) defines a local coor- 
dinate system and selects the two closest points from each sector. It guar- 
antees the neighbors to be distributed around the central point. However, 
it does not guarantee positive stencils, as the example in Fig. [5] shows. 

The stencil selected by the MPS method is shown in Fig. [6j It is minimal 
and positive, here achieving this property by selecting one point further 
away. The MPS method can be interpreted as a neighborhood criterion that 
is optimal (i.e. minimal and positive) with respect to the Laplace operator. 
A deeper discussion of various neighborhood criteria in presented in |15j . 

6 Conditions for the Existence of Positive Stencils 

We investigate when the polyhedron ([9]) is nonvoid, i.e. under which con- 
ditions positive stencils exist. We place the point of approximation in the 
origin xq = 0. A set of neighbors {x±, . . . , x m } C M. d is given, where m > k. 
In order to establish a connection between the LP space R m and the actual 
geometry space M. d we consider the dual problem, which is defined in M. k . 

Theorem 7 (Farkas' Lemma). For a real matrix A and a real vector b, 
exactly one of the following two systems has a solution: 

• A ■ x = b for some x > 0, or 

• A T ■ w > for some w satisfying b T ■ w < 0. 

3 Hybrid approaches exist that use a Delaunay mesh combined with meshfree methods. 
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Proof. The proof is given in [19] . □ 

Applying Farkas' lemma to our problem yields that system V ■ s = b has 
no solution s > 0, if and only if system V T ■ w > has a solution satisfying 
b T ■ w < 0. The i th component of V T ■ w can be written as 

(V T ■ w). = a? ■ Xi + xf • A- Xi , 
where a = (w\, . . . , Wd) T and A is the symmetric matrix 

/ \ / W 7 Wi w 5 

A = f ^ 4 ^ 3 J (2d) resp. A = i w 4 w 8 w 6 | (3d) . 

Given w (respectively a and A), we consider the quadratic form f(x) = 
a T ■ x + x T ■ A ■ x. Since A is symmetric, an orthogonal matrix S G 0{d) 
exists, such that S 7 AS = D, where D = diag (Ai, . . . , A^). In the new 
coordinates, with d = S T a, we define 

g(x) = f(Sx) = d T ■ x + x T ■ D ■ x . (11) 

If all eigenvalues Aj ^ 0, then D is reg ular. With c = -^D^d we can write 

g{x) = (x — c) T ■ D ■ (x — c) — c T ■ D ■ c , 

If one or two A, = 0, we are in a degenerate case, and stick to the represen- 
tation with d as parameter. Choosing w £ M fc arbitrarily is equivalent to 
choosing S G 0(d), A = (Ai, . . . , X d ) G M. d and c G M d (respectively d G M d ) 
arbitrarily. For any A, c G M rf define the domain 

H X c = {x G R d : > 0} . 

For a set of points X = {x±, . . . , x m } define SX = {Sx±, . . . , Sx m }. Farkas' 
lemma translates to 

Corollary 1. System V ■ s = b has no solution s > 0, if and only if 
5 G O(d), c, A G R d with £^ =1 A* < exist, such that SX C H Xc . 

In other words, no positive Laplace stencil exists, iff the set of points X 
can be transformed (via S G 0(d)), such that it is contained in the set H\ „ 

for some c, A G R d with Yfi=i \ < 0. 

Example 2. The setup in Fig. [7] shows a set of points that is completely 
contained in a domain c . Due to Cor. [H no positive stencil exists. 
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Figure 7: No positive Figure 8: Positive Figure 9: Necessary 
stencil exists stencil criterion 



Example 3. The setup in Fig. [8] shows a point setup which has a positive 
stencil solution. It is impossible to find a domain H\ _ and to rotate the 
set of points, such that all points are contained in the domain. 

We have derived a geometric condition, which is equivalent to the exis- 
tence of positive stencils. However, due to the nonlinearity in g, it is difficult 
to translate into geometric means directly. Instead, we derive a necessary 
(but not sufficient) as well as a sufficient (but not necessary) criterion on 
the point geometry for the existence of a positive Laplace stencil. To our 
knowledge the latter has not been given yet. 

6.1 A Necessary Criterion for Positive Stencils 

If V ■ s = b has a solution s > 0, then for any S G 0(d), c, A G M d with 
X^=i Aj < 0, there is a point x^ with Sxi ^ c . For the particular choice 
Ai = — 1, Xi = Vi > 1 and c\ 3> maxj \\xi\\ it follows that for any S G 0(d) 
at least one point must satisfy x\ < 0. This yields the following 

Theorem 8. If a set of points X C W^ 1 around the origin admits a positive 
Laplace stencil, then they must not lie in one and the same half space (with 
respect to an arbitrary hyperplane through the origin). 

This result is well known [3]. Due to the particular choice of A, this 
criterion is very crude, but easy to formulate in geometric means. More 
careful estimates of the condition of Cor. Q] may yield stricter criteria. 

6.2 A Sufficient Criterion for Positive Stencils 

For any c, A G W 1 with Yli=i A« < we construct a domain c D c , 
which is M. d aside from a cone centered at the origin. If for any c, A G 
S G 0(d) there is at least one point Sx^ £ c , then Sxi ^ c , thus a 
positive Laplace stencil exists. We call this criterion cone criterion. 
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Figure 10: Case Figure 11: Case Figure 12: Case Figure 13: Case 
(— ) (+") (O-)typel (0-)type2 

Theorem 9 (Cone criterion in 2d). Let c, A € M 2 with Ai + A2 < 0. There 
exists always a cone Cv defined by v ■ x > — ^==||a;||, where (3 = — 1 (a 

cone with total opening angle ^5°, the direction vector v depends on X and 
c), such that G^ c = M. d \ Cv satisfies c C G^ c . 

Proof. We show that H\ c and Cv do not intersect. Since the problem 
is invariant under interchanging coordinates, we can w.l.o.g. assume that 
A2 < 0. Including the degenerate case, three cases need to be considered: 

• Case ( ): Ai < 0, A 2 < 

The set c is the interior of an ellipse centered at c with G c . 

The vector v = — (^ci, j^cfc) is an outer normal vector. The cone Cv 
touches the ellipse only at the origin, as shown in Fig. [TU1 

• Case (H — ): Ai > 0, A 2 < 

Fig. HU shows the geometry. Define [i\ = 1^4 < 1. The domain 
H\ c is defined by g(xi,x 2 ) = ^\{x\ — 2c\X\) — [x\ — 2c 2 x 2 ) > 0. 
Due to symmetry we can assume ci,c 2 > 0. For all x £ B, where 
B = {{x\,X2)\x\ > 0, X2 < 0, \xi\ < \x2\}, the function g satisfies 

g(x) = /ii(|xi| 2 - 2ci|xi|) - (|x 2 | 2 + 2c 2 |a; 2 |) 

< (m - l)\x 2 \ 2 - 2(/iici|ari| + c 2 \x 2 \) < , 

hence c n B = 0. The domain B is a 2d cone with opening angle 
45°, where v = (±y/2-y/2, ±^2 + ^2), which proves the claim. 

• Case (0-): Ai = 0, A 2 < 

We use representation (fTT|) . Define /i 2 = |A 2 |. The domain H \ ^ is 
defined by g{x\,x 2 ) = d±x± + d 2 x 2 — (i 2 x 2 > 0. Due to symmetry we 
can assume that d\,d 2 > 0. Two subcases have to be distinguished: 
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— Case d\ ^ 0: 

The setup is shown in Fig. [T2l For all a; G B, where B = 
{{x\,X2)\x\ < 0, X2 < 0, \ xi\ < \x2\}, one has g(x\, xi) = —d\\x\\- 
d>2 \x2 1 — ^2X2 < 0, hence ^ n -B = 0. As above, the domain .B 
is a 2d cone with opening angle 45°. 

— Case d\ = 0: 

The setup is shown in Fig. [T3l The domain g(x±,X2) > is the 
set < x 2 < Any cone contained in the domain X2 < proves 
the claim. rj 

In other words, if any angle between two neighboring points (seen from 
the central point) is no more than 45°, then a positive stencil always exists. 

Remark 7. The 2d cone criterion is sharp: For any e > a point setup 
can be constructed, such that all angles are less than | + e, and a positive 
stencil does not exist. The construction is given in [15\ p. 136]. Note that 
the resulting configurations are very unbalanced. Some points are very close 
to the central point, others are far away. In practice, such extreme cases are 
typically avoided by the construction and management of the point cloud, 
yielding positive stencils also for angles significantly larger than 45°. 

Theorem 10 (Cone criterion in 3d). Let c, A G M 3 with \\ + A2 + A3 < 
0. There exists always a cone Cv defined by v ■ x > —j== ||x||, where 

P = ^/jjj(3 — y/E) ( a cone with total opening angle 33. 7° ), such that c = 

R d \ C v satisfies c cG^ c . 

Proof. The following cases need to be considered: 

• Cases ( ), (0 - -), (00-): Ai < 0, A 2 < 0, A 3 < 

As in 2d, we use representation (fTT]) . Define \ii = |Aj| Vi = 1,2,3. 
Allowing fi\ and ^2 to be zero, the domain ^ is defined by 

g(xi,x 2 ,x 3 ) = dixi + d 2 x 2 + d 3 x 3 - \i\x\ - ii 2 x\ - ^3^3 > , 

First assume that if one fa = 0, then the corresponding d{ 7^ 0. Due 
to symmetry we can w.l.o.g. assume that d\, da, d 3 > 0. For all x G B, 
where B = {(xi, X2, x 3 )\xi, X2, x 3 < 0}, the function g satisfies 

g(x 1 ,x 2 ,x 3 ) = -d^xil ~ d2\x 2 \ - d 3 \x 3 \ - /^M 2 - fJ-3\x 3 \ 2 < . 

The domain B is not a cone, but the corresponding domain from the 
case (+ H — ) is contained in it. Hence, the cone constructed in that 
case can be used here. 
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In the case that /ij = and di = 0, the geometry reduces to the 2d (or 
trivial Id) case. Since in 2d the desired estimates have been shown for 
a larger opening angle, the constructions transfer to the 3d case. 

Case (+ + -): Ai > 0, A 2 > 0, A 3 < 

Define \i\ = 1^4 < 1 and /i 2 = 1^4 < 1. Then c is defined by 

g(x 1 ,x 2 ,x 3 ) = ni(x\ - 2cxxx) + ii 2 {x\ - 2c 2 x 2 ) - (x\ - 2c 3 x 3 ) > . 

Due to symmetry we can assume c\,c 2 ,c 3 > 0. For all a; £ B, where 
B = {{xi,x 2 ,x 3 )\xi,x 2 > 0,x 3 < 0, \xi\, \x 2 \ < yj~^\x 3 \}, we have 

g(x) = /iiOi| 2 - 2ci|xi|) + n 2 (\x 2 \ 2 - 2c 2 \x 2 \) - (\x 3 \ 2 + 2c 3 |x 3 |) 
< {\{m + /U 2 ) - l)|x 3 | 2 - 2(/UiCi|xi| + [i 2 c 2 \x 2 \ + c 3 |x 3 |) < . 

Note that B is not a cone. However, a 3d cone can always be contained 
inside B. Some geometric considerations yield that the cone with max- 



imum opening angle contained inside B is given by (3 = y ^(3 — \/6) 
and v = 1 (2(V3 - y/2), 2(v / 3 - v^), l) • 

Case (+0-): A x > 0, A 2 = 0, A 3 < 

As in the preceding degenerate cases, we describe the domain by (jlll) . 
Define \x\ = |Ai| and [i 3 = |A 3 |. The case d 2 = reduces to the 2d 
case (H — ). Hence, w.l.o.g. we consider di > 0, <i 2 > 0, d 3 > 0. For all 
x £ B, where B = {(x\, x 2 , x 3 )\x±, x 2 , x 3 < 0, |xi| < \x 3 \}, we have 

g(x 1 ,x 2 ,x 3 ) = — di|ari| - d 2 \x 2 \ - d 3 \x 3 \ + ^i|xi| 2 - /x 3 |x 3 | 2 < . 

The estimate holds, since Hi < /i 3 and |a?i| < |x 3 |. As before, a 3d 
cone with desired opening angle can be contained in B. The cone from 
the case (+ H — ) can be used here. 

Case (+ - -): Ai > 0, A 2 < 0, A 3 < 

Define \i 2 = 1^4 and /i 3 = 1^4. Since \i 2 + fi 3 > 1, we assume 
w.l.o.g. jU 3 > 2- The domain _ is defined by 

g(xi,x 2 ,x 3 ) = {x\ - 2cixi) - ii 2 {x\ - 2c 2 x 2 ) - n 3 (x 2 3 - 2c 3 x 3 ) > . 
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Due to symmetry we can assume 01,02,03 > 0. For all x £ B, where 



g{x) = (|xi| 2 - 2ci\x x \) - /x 2 (H 2 + 2c 2 \x 2 \) - /i 3 (|x 3 | 2 + 2c 3 |x 3 |) 

= (|xi| 2 - fJ> 2 \x 2 \ 2 - Vz\Xz\ 2 ) -2(/XiCi|xi| + fJ. 2 C 2 \x 2 \ + C 3 |x 3 |) < 



The domain B is the same as in case (+ H — ), merely reflected at the 
xi,x% plane. Hence, a 3d cone can be placed in the same way. 

Remark 8. Unlike the 2d case, the 3d estimate is not sharp, due to the 
intermediate domain B. With significantly more algebra, it is possible to 
gain an opening angle that is a couple of degrees larger. 

Remark 9. The existence of a positive stencil implies the existence of a 
stencil. Configurations that yield an unsolvable Vandermonde system Q 
are automatically excluded by the cone criterion. 

Definition 10. We call points distributed nicely around a central point, if 
in a test cone, with opening angle given by Thm. [9] respectively Thm. [TUl 
always points are contained, for any possible direction the cone points to. 

6.3 Condition on Point Cloud Geometry 

The cone criterion guarantees positive stencils. We now provide conditions 
on the point cloud geometry and the choice of candidate points, such that 
the cone criterion is guaranteed to be satisfied. As in [10], we define 

Definition 11. Let ft C M rf be a domain and X = {x\, . . . ,x n } a point 
cloud. The mesh size h is defined as the minimal real number, such that 
& C UiLi & (xi, |) , where B (x,r) is the closed ball of radius r centered in 
x and is the closure of Q. 

We assume that a desired maximum mesh size is preserved by manage- 
ment of the point cloud, e.g. by inserting points into large holes. 

Theorem 11. Let the point cloud have mesh size h. Let 7 be the opening 
angle of the cone derived in Thm. respectively Thm. \1(K If the radius of 
considered candidate points satisfies r > sin ^/ 2 ) \ > then for every interior 
point which is sufficiently far from the boundary, a positive stencil exists. 



B = {(xi,X2,Xz)\xi > 0,X 2 ,X3 < 0, |xi|, \x 2 \ < 




x 3 |} we have 



<(!-/i 3 )NI 2 <o 
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Figure 14: Relation between neigh- Figure 15: Guaranteeing the cone 
borhood radius and mesh size criterion close to the boundary 

Proof. Having mesh size h implies that there are no holes larger in diameter 
than h, i.e. Vcc £ O 3x{ G X : \\x{ — x\\ < |. Fig. [H] shows a ball with 
radius r around the central point and a cone with opening angle 7. If the 
cone contains no point, there must be a ball of radius | which contains no 
points. The claim follows by considering the triangle (0,x,X{). □ 

The specific ratios of candidate radius to maximum hole size radius are 

r r -_ f\/4 + 2^/2 =2.61 in 2d 

h/2 > V 1 + F-j^T^ =345 in3d 

Using sharper estimates, the 3d ratio can be lowered to \/ 6 + 2\/6 = 3.30. 
In practice, point clouds are much nicer than the worst case scenario, so 
significantly smaller ratios lead to positive stencils. 

Thm. [11] is valid for any interior point which is far enough from the 
boundary, that the mesh size criterion guarantees points to lie between the 
point in consideration and the boundary. For a layer of interior points 
close to the boundary, the cone criterion can be enforced by the following 
construction (see Fig. [T5l) : First, place boundary points sufficiently dense. 
Let their maximum distance be d p . Second, ensure that every interior point 
has a minimum distance from the boundary. In 2d, this is d\> > ^d p . 

6.4 Neumann Boundary Points 

Assume the boundary dQ is C 1 around Neumann boundary points. Consider 
a local coordinate system, i.e. n = (1,0) in 2d, respectively n = (1,0,0) in 
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3d. We obtain Neumann stencils by solving the linear minimization problem 



m 

S 

mm 



V— , s.t. V- s = n, s > , (12) 



where the matrix V is given by © as 
X\ . . . x m 





( X\ .. 


2-m 


in 2d, and V = j 


yi ■ ■ 






Ui •■ 





V = ( I in 2d, and V = \ yi ... y m I in 3d. 

2/1 ■■■ ym 

We consider the 3d case. The 2d geometry is contained as a special case. 
For an easier analysis, we consider a locally convex domain, i.e. x% > Vi. 

Theorem 12. For a Neumann boundary point a positive stencil exists, iff 
the points' projections onto the normal plane do not lie all in one and the 
same half space. 

Proof. Due to Farkas' lemma, system f)12|) has no solution s > 0, iff the 
system V T -w > has a solution satisfying w x < 0, where w = (w x , w y , w z ) T . 

Let no positive stencil exist. Then w G W 1 with w x < exists, such 
that V T ■ w > 0, i.e. w x X{ + w y yi + w z Z{ > Vi This is equivalent to 

k- ( ^ % ~\ > Xj Vi, where = (j^y, jfrj)- This means that the y-z projection 

of all points lies in one and the same half space (in the direction of k). 

Conversely, assume that the y-z projection of all points lies in one and 
the same half space. Let I be the indices of all points in consideration. 
Define I p = {i £ I : Xi > 0}. Consider w.l.o.g. the case z% > \/i, where 

z t > Mi G J p . Choose w = (-1,0, ™*? ielX j ) . Then for all i G I v it holds 



w • ajj = — H > — X{ + max x% > , 

mm i6 j p ie/ 



and for all i G / \ I p one has w T ■ Xj, > 0, since Xj = 0. Thus, no positive 
stencil exists. □ 

Remark 10. Construction (|12p yields a first order accurate approximation 
of the normal derivative. Second order accuracy could be achieved, by in- 
cluding quadratic terms into V. However, in this case no positive stencil 
exists, since the condition YliLi( x i + yf) s i = cannot be satisfied. 
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Figure 16: Proper treatment Figure 17: Computational domain in 2d and 
of a crack 3d 



6.5 Treatment of Cracks 

A non-convex part of the domain f2 requires a special treatment if it is thin- 
ner than the local neighborhood radius. Fig. [16] shows a proper treatment 
of a crack. Stencils on one side of the crack must not use points on the other 
side. In the MPS method this property can be guaranteed by the following 
construction: For a central point xq E Q, only circular neighbors inside the 
star shaped core Cl Xo = {x E & : (1 — o)xq + ax E Cl Vq E [0, 1]} (bold dots 
in Fig. 1161) are considered as candidates for the linear minimization (llOp . If 
the domain is defined implicitly f2 = {x : cj>(x) < 0}, the point x does not 
lie in Q Xo if a point y E [a^O; 33 ] on the connection line satisfies (f>(y) > 0. 

7 Minimal Stencils and Matrix Connectivity 

Due to Thm. [3] and Thm. U the matrix composed of positive stencils is an 
M-matrix, if every interior and Neumann boundary point is connected to a 
Dirichlet boundary point. 

Theorem 13. Consider the Poisson problem jl| on o domain which has 
no holes (i.e. only an outer boundary). With a MPS discretization every 
interior point is connected to a boundary point. 

Proof. Assume there is a point i £ I which is not connected to a boundary 
point. Define Ij = {j E I : i is connected to j}. Every point in Jj is not 
connected to a boundary point. Hence, the set Ij C I forms an island 
inside $7 which does not reach a the boundary. Consider a point that spans 
the convex hull of Ij. It only uses points in its stencil that lie inside the 
island, hence these lie in one and the same half space, which contradicts the 
necessary condition on positive stencils given by Thm. [HJ □ 
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Figure 18: Error convergence in 2d - LSQ vs. MPS method 

Remark 11. Although Thm. [131 does not extend to interior boundaries, in 
practice the MPS works for these, given enough boundary points are placed. 

It remains to ensure that every point connects to a Dirichlet boundary 
point. Unfortunately, this cannot be concluded from the MPS method di- 
rectly. It is possible that an isolated Dirichlet point is not used in the stencils 
of nearby interior points. Note that this phenomenon can also happen on 
regular grids. A single Dirichlet point in a corner of a domain may not be 
used by regular five-point stencils. If Dirichlet data is prescribed only in 
small regions, these regions have to be equipped with a sufficient number 
of boundary points. In addition, the MPS implementation has to ensure 
that these Dirichlet points are used by nearby points. If done so, the MPS 
method guarantees to generate M-matrices. 

8 Numerical Experiments 

We investigate the numerical accuracy of the MPS method in comparison to 
a least squares approach. As test problems we consider the Poisson equation 
© in the unit box with a ball cut out Q = [0, l] d \ B{{\, . . . , |, l.l) , 0.44). 
Fig. [17] shows the computational domain in 2d and 3d. In one case, the 
boundary conditions are Dirichlet everywhere, in the other case, Neumann 
at the bottom Xd = 0, and Dirichlet everywhere else. Given g, we set 
/ = Ag and h = so ([TJ has the solution u = g. Specifically, we choose 
g(xi,x 2 ) = — (x\ sin(x 2 + 2) + x 2 sin(2xi + 1)) in 2d and g(xi,x 2 ,x 3 ) = 
7- (xi sin(:r2 + 2) + x 2 sin(2x3 + 3) + X3 sin(3xi + 1)) in 3d, with c 2 and C3 

C3 

such that max g — ming = 1. The problem is discretized by a sequence 
of point clouds. The point clouds have a uniform average density and a 
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Figure 19: Error convergence in 3d - LSQ vs. MPS method 

minimum separation [TO] of 8 = 0.05. Each point cloud is managed to satisfy 
the conditions for the existence of positive stencils, as derived in Sect. 16.31 
Since to one mesh size h, given by Def. HU many point clouds exist, we 
sample a number of experiments to obtain an average error convergence rate. 
To every point cloud we apply a weighted least squares method (Sect.HJ) and 
the MPS methods (Sect. ED, both with w(S) = 5~ 4 . 

The numerical results are shown in Fig. [18] for 2d, and Fig. [19] for 3d. 
Plotted is the error measured in the maximum norm over the mesh size h. 
Solid dots represent the all Dirichlet version of a problem, while open circles 
show the error with partial Neumann boundary conditions. The reference 
lines are of slope one and two respectively. One can observe the following: 

• Both approaches show a second order convergence rate for the pure 
Dirichlet problem, and first order convergence if Neumann boundary 
conditions are involved. While the derivation in Sect. 12.11 enforces 
only a first order accurate approximation of the Laplacian at interior 
points, point clouds tend to possess enough averaged symmetry to 
actually yield second order error convergence. On the other hand, the 
first order accurate approximation (Rem. [TO]) at Neumann boundary 
points carries through. 

• The MPS method shows a larger variation in error over the ensemble 
of experiments. One reason for this effect could be the discontinuous 
dependence on the point positions (see Rem. [I]). For two similar point 
clouds, the MPS method may select very different stencils. 

• Both methods yield roughly the same error constant. The average 
error is slightly lower with the MPS method. 
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Figure 20: CPU times for setup and solve by BiCG and AMG 
8.1 Computational Cost 

For a 3d test problem, the MPS method is compared with the LSQ method 
in terms of computational cost. For a sequence of point clouds, Fig.l20lshows 
the CPU times for the setup of the system matrix (left plot), the solution of 
the arising system with a BiCG scheme (center plot), and the solution with 
an algebraic multigrid (AMG) scheme (right plot). The latter is performed 
using SAMG |14j by the Fraunhofer Institute for Algorithms and Scientific 
Computing. 

As expected (Sect. I5.ip . the cost of setting up the system matrix is 
roughly equal for MPS and LSQ method. In fact, MPS is slightly slower 
with the used simplex method. However, more efficient linear programming 
methods may turn the tide towards the MPS method. On the other hand, 
the cost for solving the large linear system is significantly reduced by the 
MPS method. The speedup factor equals the factor in sparsity, i.e. the MPS 
approximation does not modify the convergence rate. While the AMG solver 
shows a cost roughly linear in the number of unknowns, solvers further away 
from optimal effort (like BiCG) will greatly benefit from the sparsity of the 
MPS approximation as the number of unknowns increases. 

9 Conclusions and Outlook 

We have presented a meshfree approach that constructs minimal positive 
stencils for the Laplace operator on a cloud of points. We have shown that 
under moderate assumptions on the local resolution of the point cloud, posi- 
tive stencils always exist. The method approximates the Poisson equation by 
M-matrices, which are optimally sparse. Both properties are beneficial, and 
they are in general not met by classical least squares approaches. Numer- 
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ical tests show that with the presented approach, the approximation error 
roughly equals the error of classical approaches. The computational cost to 
construct the new approximation is comparable to the cost of least squares 
methods. On the other hand, for solving the arising linear system, both 
cost and memory requirements are reduced significantly due to the optimal 
sparsity. 

An efficient solution of the linear programs is a crucial point in the 
presented method, worth a deeper analysis. The application to particle 
methods shall be investigated. While the presented approach can stand as 
a method of its own, it may also increase the efficiency of other approaches. 
The minimal stencils can be augmented by additional neighbors and the 
final stencil be computed by a least squares method. For instance, a local 
neighborhood radius can be based on the farthest minimal stencil point, thus 
increasing sparsity and adding local adaptivity to existing meshfree codes. 
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